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■ Abstract. In this paper, Lipschitz univariate constrained global optimization problems where 
both the objective function and constraints can be multiextremal are considered. The con- 
strained problem is reduced to a discontinuous unconstrained problem by the index scheme 
without introducing additional parameters or variables. A Branch-and-Bound method that 
does not use derivatives for solving the reduced problem is proposed. The method either 
determines the infeasibility of the original problem or finds lower and upper bounds for the 
global solution. Not all the constraints are evaluated during every iteration of the algorithm, 
providing a significant acceleration of the search. Convergence conditions of the new method 
are established. Test problems and extensive numerical experiments are presented. 
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1. Introduction 



H 

■ Global optimization problems arise in many real-life applications and were 

intensively studied during last decades (see, for example, (Archetti and Scho- 
en, 1984; Bomze et al., 1997; Breiman and Cutler, 1993; Evtushenko, 1992; 
Floudas and Pardalos, 1996; Horst and Pardalos, 1995; Horst and Tuy, 1996; 
Locatelli and Schoen, 1999; Lucidi, 1994; Mladineo, 1992; Pardalos and 
Rosen, 1990; Pinter, 1996; Strongin, 1978; Sun and Li, 1999; Torn and Zilin- 
skas, 1989; Zhigljavsky, 1991), etc.). Particularly, univariate problems attract 
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attention of many authors (see (Calvin and Zilinskas, 1999; Hansen and Jau- 
mard, 1995; Lamar, 1999; Locatelli and Schoen, 1995; MacLagan, Sturge, 
and Baritompa, 1996; Pijavskii, 1972; Sergeyev, 1998; Strongin, 1978; Wang 
and Chang, 1996)) at least for two reasons. First, there exist a large number of 
applications where it is necessary to solve such problems (see (Brooks, 1958; 
Hansen and Jaumard, 1995; Patwardhan, 1987; Ralston, 1985; Sergeyev et 
al., 1999; Strongin, 1978)). Second, there exist numerous schemes (see, for 
example, (Floudas and Pardalos, 1996; Horst and Pardalos, 1995; Horst and 
Tuy, 1996; Mladineo, 1992; Pardalos and Rosen, 1990; Pinter, 1996; Stron- 
gin, 1978)) enabling to generalize to the multidimensional case the mathe- 
matical approaches developed to solve univariate problems. 
In this paper we consider the global optimization problem 

mm{f{x) :xe[a,b], gj{x) <0, l<j<m}, (1) 

where f{x) and gj{x), 1 < j < m, are multiextremal Lipschitz functions (to 
unify the description process we shall use the designation gm+i{x) = f(x)). 
Hereinafter we use the terminology "multiextremal constraint" to highlight 
the fact that the constraints are described by multiextremal functions gj (x) , 1 < 
j < m, in the form (1) (of course, the same subregions of the interval [a,b\ 
may be defined in another way). In many practical problems the order of 
the constraints is fixed and not all the constraints are defined over the whole 
search region [a,b\ (if the order of the constraints is not a priori given, the user 
fixes his/her own ordering in a way). In the general case, a constraint gj+\ {x) 
is defined only at subregions where gj (x) < 0. We designate subdomains of 
the interval [a,b] corresponding to the set of constraints from (1) as 

Qi = [a,bl Qj+i = {x(EQj:gj{x)<0}, 1 < ; < m, (2) 
We introduce the number M such that 

GmT^O, QM+l=QM+2--- = Qm+l=^- (3) 

If the feasible region of the problem (1) is not empty then Qm+\ 7^ and 
M = m + 1 . In the opposite case M indicates the last subset Qj from (2) such 
that Qj / 0. 

We suppose in this paper that the functions gj{x),\ < j <m + l, satisfy 
the Lipschitz condition in the form 

\gj{x')-gj{x")\<Lj\x'-x"\, x',x"eQj, l<y<m+l. (4) 

where the constants 

0<L^<oo, l<j<m + l, (5) 
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are known (this supposition is classical in global optimization (see (Hansen 
and Jaumard, 1995; Horst and Tuy, 1996; Pijavskii, 1972)), the problem of 
estimating the values Lj,l < j < m + 1, is not discussed in this paper). Since 
the functions gj{x),l < j < m, are supposed to be multiextremal, the sub- 
domains Qj,2 < j <m + l, can have a few disjoint subregions each. In the 
following we shall suppose that all the sets Qj,2 < j <m+l, either are empty 
or consist of a finite number of disjoint intervals of a finite positive length. 

In the example shown in Fig. 1, a) the problem (1) has two constraints 
g\{x) and g2{x). The corresponding sets Qi = [a,b],Q2, and Qt, are shown. It 
can be seen that the subdomain Q2 has three disjoint subregions and the con- 
straint giix) is not defined over the subinterval [c,d]. The objective function 
f{x) is defined only over the set Q3. 

The problem (1) may be restated using the index scheme proposed origi- 
nally in (Strongin, 1984) (see also (Strongin and Markin, 1986; Strongin and 
Sergey ev, 2000)). The index scheme does not introduce additional variables 
and/or parameters by opposition to classical approaches in (Bertsekas, 1996; 
Bertsekas, 1999; Horst and Pardalos, 1995; Horst and Tuy, 1996; Nocedal 
and Wright, 1999). It considers constraints one at a time at every point where 
it has been decided to calculate gm+i{x). Each constraint gi{x) is evaluated 
only if all the inequalities 

gj{x) < 0, 1 < < /, 

have been satisfied. 

In its turn the objective function g^+i {x) is computed only for that points 
where all the constraints have been satisfied. 

Let us present the index scheme. Using the designations (2), (3) we can 
rewrite the problem (1) as the problem of finding a point xl^ and the corre- 
sponding value glf such that 

8*M = 8m{x*m) = mm{gM{x) : x £ Qm}. (6) 

The values coincide with the global solution of the problem (1) if 

M = m + 1, i.e. when the original problem is feasible. We associate with every 
point of the interval [a, b] the index 

v = v{x), l<v<M, 

which is defined by the conditions 

gj{x)<o, i<y<v-i, gvW>o, (7) 

where for v = m -I- 1 the last inequality is omitted. We shall call trial the 
operation of evaluation of the functions gj{x), 1 < j < v(x), at a point x. Let 
us introduce now an auxiliary function (p(x) defined over the interval [a,b] as 
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follows 

/ \ / \ f 0) if v(x) < m + 1 

cpW=gvwW-|^^^^^ ifv(x)=m+l 

where is the solution to the problem (1) and to the problem (6) in 
the case M = m+\. Due to (6), (8), the function (p(x) has the following 
properties: 

i. (p(x) > 0, when v(x) < m + 1; 

ii- > 0, when v(x) = m + 1; 
iii. (p(;c) = 0, when v(x) = m + 1 and g^+i (x) = g^^j. 

In this way the global minimizer of the original constrained problem (1) 
coincides with the solution x* of the following unconstrained discontinuous 
problem 

(p(x*) = raia{(^{x) : x G [a,b\}, (9) 

in the case M = m + \ and gm+\{x*) = g*„^i- Obviously, the value g^^j 
used in the construction (8) is not known. Fig. 1 b) shows the function (p(x) 
constructed for the original problem from Fig. 1 a). 

Numerical methods belonging to the class of information algorithms based 
on probabilistic ideas have been proposed for solving the problem (9) in 
(Strongin, 1984; Sergeyev and Markin, 1995; Strongin and Markin, 1986; 
Strongin and Sergeyev, 2000). 

In this paper a new method called Index Branch-and-Bound Algorithm 
(IBBA) is introduced for solving the discontinuous problem (9). The next sec- 
tion shows that, in spite of the presence of unknown points of discontinuity, it 
is possible to construct adaptively improved auxiliary functions (called by the 
authors index support functions) for the function (p(x) and to obtain lower and 
upper bounds for the global minimum. The computational scheme of the new 
method is described in Section 3. Convergence conditions of the algorithm 
are established in Section 4. Section 5 contains wide computational results 
showing quite a promising behaviour of the new algorithm. Finally, Section 
6 concludes the paper. 



2. Discontinuous index support functions 

It has been shown in (Pijavskii, 1972) that lower and upper bounds can be 
found for the global solution F* of the problem 

F* =mm{F{x):xe[a,b]}, (10) 

where 

\F{x')-F{x")\<Lf\x' -x" \, x',x"e[a,b], (II) 
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through sequential updating of a piece-wise linear support function 

\|/(x)<F(x), x^[a,b], (12) 

if the Lipschitz constant < L/r < oo is known. The algorithm proposed in (Pi- 
javskii, 1972) improves the support function during every iteration by adding 
a new point where the objective function F{x) is evaluated. This procedure 

enables to draw the support function closer to the objective and, therefore, 
to decrease the gap between the lower and upper bounds. Let us show that 
by using index approach it is possible to propose a procedure allowing to 
obtain lower and upper bounds for the solution g*„^i. In order to induce the 
exhaustiveness of the partitioning scheme in the further consideration it is 
supposed that constants Kj such that 

Lj<Kj<oo, \<j<m + \, (13) 

are known. The case Lj = Kj is discarded from the further consideration 
because in the algorithm of Pijavskii it leads to a possibility of generation 
of a new point coinciding with one of the points previously generated by the 
method. 

Suppose that k trials have been executed at some points 

a=xo <xi < ... <Xi < . .. <Xk = b (14) 

and the indexes v,- = v(x;),0 < / < ^, have been calculated in accordance with 
(7). Since the value gl^_^_^ from (8) is not known, it is not possible to evaluate 
the function cp(;<:) for the points having the index m + 1. In order to overcome 
this difficulty, we introduce the function (p^(x) which is evaluated at the points 
Xi and gives us the values Zi = (Pk{xi),0 <i<k, as follows 

/ \ / \ f ifv(x) <m+l ,,,, 

9.W=gv(x)W-|2. ifv(x)=m+l ^^^^ 

where the value 

= min{g^+i (x/) : < i < A:, V,- = m + 1}. (16) 

estimates from (8). It can be seen from (8), (15), and (16) that (Pk{xi) = 
(p(x;) for all points x/ having indexes v(x;) < m + 1 and 

< < 

if v(x/) = m + 1 . In addition, 

cp-tW < 0, x€{x: gm+i (x) < Zl}. (17) 
During every iteration the trial points Xi,0 <i <k, form subintervals 
[xi-i,Xi] (Z[a,b\, 1 < j < k, 
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and every point x; has its own index v, =v{xi),0 <i <k, calculated in accor- 
dance with (7). Then, there exist the following three types of subintervals: 

i. intervals such that v,_i =v,; 

ii. intervals such that v,_i < V,; 

iii. intervals such that V,_i >V,-. 

The bounding procedure presented below constructs over each interval 
[;c,_i,x,] for the function (pk{x) from (15) a discontinuos index support func- 
tion "^i{x) with the following properties 

(^) <^k{x), xe [xi- 1 , X,-] n 

where 

v7 = max{v(x,_i),v(x;)}. 

Note that the introduced notion is weaker than the usual definition of a support 
function (cf. (12)). In fact, nothing is required with regard to behaviour of 
\|/;(x) over [x,_i,x,] \2v7 and \|/i(x) can be greater than (p;t(x) on this subdo- 
main. 

Let us consider one after another the possibilities (i)-(iii). The first case, 
v,_i = Vj, is the simplest one. Since the indexes of the points x;_i ,x; coincide, 
the index support function is similar to that one proposed in (Pijavskii, 1972). 
In this case, due to (4), (13), and (Pijavskii, 1972), we can construct for (pit(x) 
the index support function \|/,(x),x G [x,_i,x;], such that 

> W> ^ ^ [xi-\,Xi] n gvi, 

where the function \|//(x) (see (15), (16)) has the form 

\|//(x) = max{^v,(-^i-i) --Kv,- I Xi-i -X \,gyXxi)-K^, | x,- -x |} (18) 

in the case v,_i = v, < m + 1 and the form 

\|/,(x) = max{g^+i (x;_i ) -Zl- K^+i | X;_l - x | , 

8m+i {xi) -zl- K„t+i I X, - X I } (19) 

in the case V/_i = v,- = m + 1; the constants K^,. are from (13). In both cases 
the global minimum Rj of the function \|/, (x) over the interval [x/_i ,x,] is 

Ri = 0.5(z/_i +zi - Ky,i^i-Xi-i)), (20) 

and is reached at the point 

ji = 0.5(x,-_i +Xi - {zi - Zi-i)/Ky,). (21) 
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This case is illustrated in Fig. 2 where the points Xi,Xi+i, ends of the interval 
have the indexes v,- = V/+i = j + l <m+l.In this example 

gi (xi) <0,...,gj (xi) < 0, g j+i (xi) > 0, 

gi{xi+\) <0,...,gj{xi+i) <0, gj+i{xi+i) >0, 

Zi = (Pk{Xi) =gj+l{Xi), Zi+l = (?k{Xi+\) =gj+\{Xi+\)- 

The values and are also shown. The interval [x,_2,x/_i] in the same 
Figure illustrates the case v,_2 = v,_i = j. 

The second case is v,_i < v,. Due to the index scheme, this means that 
the function (pk{x) has at least one point of discontinuity (O over the interval 
(see an example in Fig. 2) and consists of parts having different 
indexes. To solve the problem (9) we are interested in finding the subregion 
having the maximal index M from (3). The point x; has the index v, > v,_i 
and, due to (15), we need an estimate of the minimal value of the function 
<Pk{x) only over the domain [;iC;_i,;c,] n Qy-. The right margin of this domain 
is the point x; because it is the right end of the interval [x,_i ,x,] and its index 
is equal to v,. It could be possible to take the point as an estimate of the 
left margin of the domain ,x,] fl Qy. but a more accurate estimate can be 
obtained. 

It follows from the inequality v,_i < V,- that 

Zi-\ = (?k{xi-i) = gv,_t (xi-i) > 0, ^v,_i (xi) < 0. 
The function gy._^ (x) satisfies the Lipschitz condition, thus 

gv,_i (x) > 0, xe [xi-uyY ) 1^ 2v,_i , 
where the point is obtained from (18) 

yy = x,_ 1 + Zi- 1 / ^v,_ 1 • (22) 

An illustration of this situation is given in Fig. 2 where the point CO G [y,",^^!] 
is such that gvi_i (k>) = and 

[xi-i,Xi\ n Qy._^ = [xi-i, co] , [xi-i,Xi] n Qy. = [co,x,-] , 

[x,-_i,yr]n(2v,._i = [xi-i,y^]. 

Therefore, the function gyj{x) can be defined at most over the interval 

[y,~,x,] and the point yy can be used as an estimate of the left margin of the 
set [x,_i,x,] n Qy. for finding a lower bound for the function (pyt(x) over this 
domain. The corresponding index support function \|/, (x) in this case has the 
form 

\|/,-(x) =z,--ii:v; |x;-x| (23) 
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and, therefore, 

min{\|/i(x) : x € [yj ,Xi\} < min{\|/,(x) : x € [y'i',Xi] n Qv,}- 
This minimum is located at the point y'^ and can be evaluated as 

Ri = Zi - Ky. {xi -yY)= Zi - Ky. (x,- - i - Zi- 1 / Ky._^ ) . (24) 

Let us consider the last case v,_i > v,- being similar to the previous one. 

The point x; has the index v,- < v,_ i and, due to the index scheme, we need an 
estimate of the minimal value of the function (pjt(x) over the domain [x,_i ,x,] fl 

ev,_.. 

Since we have v,_i > v,-, it follows 

Zi = (Pi:(x;) = gyXxi) > 0, gv,(-^'-l) < 0- 

The function gy. (x) satisfies the Lipschitz condition and, therefore, 
gv,W>0, xe{yt,xu]f^Qy., 

where 

yj =Xi-Zi/Ky,. (25) 

Thus, the function gy._^ (x) can be defined at most over the interval [x;_i , j^]. 
The corresponding index support function 

\|/, (x) = Zi- 1 - ifvi_i \xi-i-x\. (26) 

It is evident that 

min{\|/,(x) :xG [x;_i,yt]} < min{\|/,(x) : x G [xi-l,yf]f^Qy^_^}. 
It is reached at the point yf and can be calculated as 

Ri=Zi-i -Ky^_^{yl -Xi^l) =Zi-l -Ky__^{Xi-Xi^l-Zi/Ky^). {11) 

This case is illustrated in Fig. 3. The points x;_i,x,- have the indexes V;_i = 
J + 2 < m + 1 , V; = j. This means that 

g\{Xi-l) <0,...,gj+l{Xi-l) <0, g;+2(Xi__l) >0, 

gi {xi) < 0, . . . ,gj-i (x,) < 0, gj{xi) > 0. 
The values Zi-i and Zi are evaluated as follows 

Zi-l = (?k{Xi-l) = gj+2{Xi-l), Zi = (?k{Xi) = gj{Xi). 

Fig. 3 presents a more complex situation in comparison with Fig. 2. In 
fact, 

[xi-uyt]<^Qvi-i = [xi-i,(o\\{Qj+i n [xi-i,xi]}. 
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The existence of the subregion Qj+i n [x,_i,;c,] cannot be discovered by 
the introduced procedure in the current situation because only the information 

Xi- 1 , V/_ 1 , Kyi_ ^ , Zi- 1 , Xi,Vi, Km. , Zi 

regarding the function (^k{x) over is available. This fact is not rel- 

evant because we are looking for subregions with the maximal index M, 
i.e. subregions where the index is equal to j+ \ are not of interest because 
M > + 2 since v,_i = j + l. 

Now we have completed construction of the function In all three 

cases, (i) - (iii), the value /?; being the global minimum of over the 
interval has been found (hereinafter we call the value Ri character- 

istic of the interval [x, i,x;]). It is calculated by using one of the formulae 
(20),(24), or (27) and is reached at the points from (21), yj is from (22), or 
yf from (25), correspondingly. 

If for an interval [x,_i ,x,] a value Ri > has been obtained then, due to the 
index scheme, it can be concluded that the global solution x*^_^-^ ^ 
For example, in Fig. 2 the intervals 

have positive characteristics and, therefore, do not contain the global mini- 
mizer. 

Let us now consider an interval [x,_i,Xi] of the type (iii) having a negative 
characteristic (see, for example, the interval [x,_i,x;] from Fig. 3). The value 
Ri < has been evaluated at the point yf as the minimum of the function 
\|/,(x) from (26). Since Zi-\ = \|/,(x,_i) > and = "^i{yf) < 0, a point 
X G [x,_i,j;+] such that \J/i(x) = can be found. It follows from (22) that 
X = yi • Thus, the subinterval [y,",)'/'] is the only set over [x,_i,x,] where the 
function cpi:(x) can be less than zero and where, therefore, the global solution 
x*^^i can possibly be located. By analogy, it can be shown that when < 
in the cases (i) and (ii), the interval [j,^,}'/^] is again the only subinterval of 
[x,_i,x,] where the global solution can possibly be located. 

The Index Branch-and-Bound Algorithm (IBBA) proposed in the next sec- 
tion at every (/: + l)th iteration on the basis of information obtained during the 
previous k trials constructs the function (Pfe(x) and the index support functions 

{x),'^ <i< k. Among all the intervals [x,_ i , x,] , 1 < / < it finds an interval 
t with the minimal characteristic R(, and chooses the new trial point x^+^ 
within this interval as follows 

f 0.5{yi+y^), Vt-i =Vt 
x^+' = l 0.5{yr+x,), Vt-i<Vt (28) 
^ 0.5(x,_i+y+), V,_i>V, 

Note that for intervals having V(_i = the new trial point x^+^ coincides 
with the Pijavskii point y,, i = t, from (21). Thus, the new algorithm at every 
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iteration updates the function (pyt(^) making it closer to cp(x) trying to improve 
the estimate of the global minimum 

3. Description of the algorithm 

Let us describe the decision rules of the IBBA. The algorithm starts with 
two initial trials at the points = a and = b. Suppose now that: a search 
accuracy 8 has been chosen; k trials have been akeady done at some points 
x^, . . . ,;c*; their indexes and the value 

= max{v(x') : < ? < A:} (29) 

have been calculated. Here the value estimates the maximal index M from 
(3). 

The choice of the point x^"*" ^ ,k> 1 , at the (A: + 1 )-th iteration is determined 
by the rules presented below. 

Step 1. The points x*^, . . . of the previous k iterations are renumbered by sub- 
scripts in order to form the sequence (14). Thus, two numerations are 
used during the work of the algorithm. The record x^ means that this point 
has been generated during the A:-th iteration of the IBBA. The record Xk 
indicates the place of the point in the row (14). Of course, the second 
enumeration is changed during every iteration. 

Step 2. Recalculate the estimate from (16) and associate with the points x, the 
values Zi = (pjt(x,),0 <i <k, where the values (p,t(x,) are from (15). 

Step 3. For each interval [x;_i,x,], 1 < j < A:, calculate the characteristic of the 
interval 

{0.5 (Zi_ i+Zi- {Xi - X;_ 1 ) ) , V;_ 1 = V; 

Zi-Ky_{xi-Xi-i -Zi-i/Ky^^,), Vi-i < V; (30) 
Zi-\ - K^i_i {Xi-Xi-i - Zi/Ky,.), V,_i > V; 

Step 4. Find the interval number t such that 

t = min{argmin{/?,- : 1 < ? < k}}. (31) 

Step 5. (Stopping Rule) If Rf > 0, then Stop (the feasible region is empty). Oth- 
erwise, if 

xt-xt-i>e (32) 

go to Step 6 (8 is a preset accuracy and t is from (31)). In the opposite 
case. Stop (the required accuracy has been reached). 
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Step 6. Execute the {k + l)-th trial at the point from (28), evaluate its index 
v(x*+^) and the estimate M*^+^ and go to Step 1. 

In the following section we will gain more insight the method by estab- 
lishing and discussing its convergence conditions. 



In this section we demonstrate that the infinite trial sequence {x*^} generated 
by the algorithm IBBA (8 = in the stopping rule) converges to the global 
solution of the unconstrained problem (9) and, as consequence, to the global 
solution of the initial constrained problem (1) if it is feasible. In the opposite 
case the method establishes infeasibility of the problem (1) in a finite number 
of iterations. 

In Lemma 1, we prove the exhaustiveness of the branching scheme. The 
convergence results of the proposed method can be derived as a particular 
case of general convergence studies given in (Horst and Tuy, 1996; Pinter, 
1996; Sergeyev, 1999). We present a detailed and independent proof of these 
results in Theorems 1 and 2. 

LEMMA I. Let x be a limit point of the sequence {x^} generated by the 
IBBA with £ = in the stopping rule (32), and let i = i{k) be the number of 
an interval [x,(;t)-i)^((it)] containing this point during the k-th iteration. Then 



Proof: During the current k-t\i iteration an interval [x(_i,X(] is chosen for 
subdivision. Due to the decision rules of the IBBA and (22), (25), this means 
that its characteristic Rt <0 and the point x*^"*"^ from (28) falling into the 
interval {xt-\,Xt) can be rewritten as follows 



4. Convergence conditions 




(33) 




Vj-l =v, 

Vf-i <y, 

Vf-i >Vf 



(34) 




(35) 



Let us show that the following contracting estimate 



max{x; — x' 



.k+\ 



0.5(1+ max{Lv,_i / Ky,_^ ,LyjKy,}){x, -x,-i) 



(36) 
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holds for the intervals (35), where 

0.5 < 0.5(1 +m^x{L,,_jK^,_„L^jKy,}) < 1. (37) 

Let us consider three cases. 

i. In the first case \t-i = v?, and Zt-\ = gy,{^t-\), Zt = gv,{xt)- It follows 
from (11), (13) that 

I Zt-Zt-1 |<Lv,(x,-x,_i) 

and, due to (21), we have 

max{x; - /+\/+^ - 1 } < 0.5 ( 1 + Lv,_i /Ky,_^ ){xt-Xt-i). 

Thus, (36) and (37) have been established. 

ii. In the second case, Vf_i < V(, and, therefore, due to the index scheme, 
Zt-i = gv,_i (xt-i) > and gv,_i (xt) < 0. From this estimate and the obvious 
relation 

8\,-i{xt) > Zt-i -Lv,_i(x, -x,-i) 

we obtain 

Zf_i -Lv,_|(x, < 0. (38) 
Since Zt-i > 0, it follows from (38), (28), and (13) that 

-xt-i = 0.5{xt -xt-i +Zt-i/Ky,_i) < 

0.5{xt -xt-i +l^,_jK^,_^{xt -xt-i)). (39) 
Let us now estimate the difference Xt — ;c*+^ 

X, = 0.5(x, -Xt_r-z,-i/K^,_,) < 0.5(x, -x,_^). (40) 

Obviously, the estimate (36) is the result of (39) and (40). 

iii. The case Vj i > is considered by a complete analogy to the case (ii) 
and leads to estimates 

-xt-x <0.5(x,-x,_i), 

Xt -x*"*"^ < 0.5(x, -x,_i +LyjK^^{xt -x,_i)). 

To prove (37) it is enough to mention that Ly^_^,Ky^_^,Ly,, and ^Ty, are 
constants and (13) takes place for them. The result (33) is a straightforward 
consequence of the decision rules of the IBBA and the estimates (36), (37). ■ 

THEOREM 1 . If the original problem (1) is infeasible then the algorithm 
stops in a finite number of iterations. 
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Proof: If the original problem (1), (4) is infeasible then the maximal index M 
over the interval [a,b] is less than m + 1. In this case (see (7), (8), and (15)) 

(p(x) = (pi.(x) > 0, xe[a,b]. 

On one hand, due to (4), (13), the linear pieces of the index support functions 
^i{x),i <i<k, from (18), (23), and (26) constructed by the algorithm have a 
finite slope. On the other hand, Lemma 1 shows that the length of any interval 
containing any limit point goes to zero. 

Thus, it follows from our supposition regarding the sets Qj,l < j <m + l, 
being either empty or consisting of a finite number of disjoint intervals of a 
finite positive length and the formulae (37), (30), and (31) that there exists a 
finite iteration number such that a characteristic Rt(N) > will be obtained 
and the algorithm will stop. ■ 

Let us now consider the case when the original problem (1) is feasible. 
This means that M = m + 1 in (6), (8). Let us denote by X* the set of the 
global minimizers of the problem (1) and by X' the set of limit points of the 
sequence {x*} generated by the IBBA with e = in the stopping rule (32). 

THEOREM 2. If the problem (1 ) is feasible then X* = X'. 

Proof: Since the problem (1) is feasible, the sets Qj, 1 < j < m + 1, are not 
empty and therefore, due to our hypotheses, they consist of a finite number 
of disjoint intervals of a finite positive length. This fact together with 8 from 
(32) equal to zero leads to existence of an iteration q during which a point 

having the index m + 1 will be generated. Thus, (see (15), (16)) the first 
value z; = corresponding to the point )fl will be obtained and during all 
the iterations k> q there will exist at least two intervals having negative 
characteristics (see (30)). 

Let us return to the interval [x;_i,x,] from Lemma 1 containing a limit 
point X G X'. Since it contains the limit point and the trial points are chosen 
by the rule (31), its characteristic should be negative too for all iterations 
k> q. Then, by taking into consideration the facts that Zi > 0, 1 < j < A:, (see 
(15)) it follows from Lemma 1 and (30), (31) that 

Um%)=0. (41) 

We can conclude from (41) and < 0,/: > ^, that 

lim (p^(x) = (p(x) = 0. (42) 

Let us consider an interval [xj{j^-^_i,Xj(j^)\ containing a global minimizer 
X* G X* during an iteration k> q. Ai first, we show that there will exist an 
iteration number c>q such that ^j(c)-\ =m + \or Vy(c) = m + 1. If trials will 



JOGO_indexBB.tex; 15/01/2013; 13:41; p. 13 



14 



Ya. D. Sergeyev, D. Famularo, and P. Pugliese 



fall within the interval due to the decision rules of the IBBA, 

such a trial will be generated. Suppose that trials will not fall into this interval 
and 

r = max{v (xji^k) -i),^{xj{k))} On+l. 
The point x* is feasible, this means that 

X* e [a,p] = [xj(k)-i,Xj(k)]<^Qm+h 
where the interval [a, P] has a finite positive length and 

gl{x)<0, \<l<m, xG[a,P]. 
We obtain from these inequalities that, due to (30) and (13), the characteristic 

Rj(^k) < min{gr (^) : ^ G [a, P] } < 0. (43) 

Since trials do not fall at the interval [xy(yt)-i,-^;(/t)]» it follows from (30) that 
is not changed from iteration to iteration. On the other hand, the charac- 
teristic when k ^ °°. This means that at an iteration number k' > q 
the characteristic of the interval , / = i{k'), will not be minimal. Thus, 
a trial will fall into the interval [xj-i,Xj]. The obtained contradiction proves 
generation of a point 

x''G[a,P], c>k', v{j^)=m + l. 

We can now estimate the characteristic Rj(^f,^ of the interval 
containing the global minimizer x* G X* during an iteration k> c. We have 
shown that at least one of the points Xj(^)_i,Xj(y(.) will have the index m + l. 
Again, three cases can be considered. 

In the case Vy_i = Vy = m + 1 it follows from (4), (15) that 

Zj-l - < Ln+\{x* -Xj-l). 

From (17) we have (p;t(x*) < and, therefore, 

Zj-\ <Lm+l{x*-Xj-\), 

Analogously, for the value Zj it follows 

Zj < Lm+l{Xj —X*). 

From these two estimates we obtain 

Zj + Zj-\ <Lm+i{xj-xj^i). 
By using (13), (20), and the last inequality we deduce 

Rj{k) < {Lm+l-Km+l){xj(k) -Xj(^k)-l) < 0- (44) 



JOGO_indexBB.tex; 15/01/2013; 13:41; p. 14 



Index Branch-and-Bound Algorithm for Global Optimization with Multiextremal Constraints 15 



Analogously, it can be seen from (22) and (24) that in the case Vj-i < Vj the 
estimate 

Rj{k) < {Lm+i - K„+i) {xj^k) - yj^j^-^ ) < (45) 
takes place because 

Zj < Lm+l{Xj -X*) < L„+l{Xjlj,) -y7(fc))- 

For the case > Vj (see (25) and (27)) we have 

Rm < (Lm+i - Km+i){y^^k) -""m-i) < 0- (46) 

It follows from (43) - (46) that the characteristic of the interval [;cj_i,A:y] 
containing the global minimizer x* will be always negative. Assume now, 
that X* is not a limit point of the sequence {x'^}, then there exists a number P 
such that for a\\k>P the interval , j = j{k), is not changed, i.e. new 

points will not fall into this interval and, as a consequence, its characteristic 
Rj^k) will not change too. 

Consider again the interval from Lermna 1 containing a limit 

point X € X'. It follows from (41) and the fact that is a negative constant 
that there exists an iteration number A^^ such that 

Rm<Ri{N)- 

Due to decision rules of the IBBA, this means that a trial will fall into the 
interval [xj-i,Xj]. But this fact contradicts our assumption that x* is not a 
limit point. 

Suppose now that there exists a limit point x EX' such that x ^X*. This 
means that (p(jc) > (p(x*),x* G X*. Impossibility of this fact comes from (41), 
(42), and the fact of x* € X'. ■ 

We can conclude that if the algorithm has stopped and has not established 
that Qm+i = then the following situations are possible: 

i. If < m + 1, then this means that the accuracy 8 was not sufficient for 
establishing the feasibility of the problem; 

ii. If = m + 1 and all the intervals [xp-i,Xp] such that 

max{Vp_ 1 , Vp } < m + 1 (47) 

have positive characteristics then, we can conclude that the global mini- 
mum z* of the original problem (1) can be bounded as follows 

where the value is from (16) and is the characteristic correspond- 
ing to the interval number t = t{k) from (31). 
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iii. If M* = m + 1 and there exists an interval [xp-i,Xp] such that Rp <0 
and (47) takes place then, the value can be taken as an upper bound 
of the global minimum z* of the original problem (1). A rouge lower 
bound can be calculated easily by taking the trial points xt such that 
v(x,) = m + 1 and constructing for f{x) the support function of the type 
(Pijavskii, 1972) using only these points. The global minimum of this 
support function over the search region [a, b] will be a lower bound for z* . 
A more precise lower bound can be obtained by minimizing this support 
function over the set 

U[x,_i,x,], Ri<0, l<i<k. 

We do not discuss here the pecuharities of the implementation of the 
IBBA. Let us make only two remarks. First, it is not necessary to re-calculate 
all the characteristics during Step 3 but it is sufficient to do this operation 
only for two new intervals generated during the previous iteration. Second, 
as it follows from the proofs of Theorems 1 , 2, it is possible to exclude from 
consideration all the intervals having positive characteristics. 

5. Numerical comparison 

The IBBA algorithm has been numerically compared to the method (indicated 
hereinafter as PEN) proposed by Pijavskii (see (Pijavskii, 1972; Hansen and 
Jaumard, 1995)) combined with a penalty function. The PEN has been chosen 
for comparison because it uses the same information about the problem as the 
IBBA - the Lipschitz constants for the objective function and constraints. 

Ten differentiable and ten non-differentiable test problems introduced in 
(Famularo, Sergeyev, and Pugliese, 2001) have been used. In addition, the 
IBBA has been applied to one differentiable and one non-differentiable in- 
feasible test problem from (Famularo, Sergeyev, and Pugliese, 2001). Since 
the order of constraints can influence speed of the IBBA significantly, it has 
been chosen the same as in (Famularo, Sergeyev, and Pugliese, 2001), without 
determining the best order for the IBBA. The same accuracy 8 = 10~^ {b — a) 
(where b and a are from (1)) has been used in all the experiments for both 
methods. 

In Table I (Differentiable problems) and Table II (Non-Differentiable prob- 
lems) the results obtained by the IBBA have been summarized and the columns 
in the Tables have the following meaning: 

- the columns XIBBA and FIBBA represent the estimate to the global 
solution {x*,f{x*)) found by the IBBA for each problem; 

- the columns Ng^, Ng^, Ng^ represent the number of trials where the con- 
straint gi,l < j < 3, was the last evaluated constraint; 
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Table I. Results of the experiments executed by the IBBA with the differentiable problems. 



Problem 


XIBBA 


FIBBA 








Nf 


Iterations 


Eval. 


1 




— / .D 12/0549 








1 o 
IJ 


Zi 


Jo 


z 






zUo 






zi 


Lit 


z4s 


3 


-5.99182849 


-2.94266082 


40 






22 


62 


84 


4 


2.45965829 


2.84080900 


622 


156 




175 


953 


1459 


5 


9.28501542 


-1.27484676 


8 


14 




122 


144 


402 


6 


2.32396546 


-1.68515824 


14 


80 




18 


112 


228 


7 


-0.77473979 


-0.33007410 


35 


18 




241 


294 


794 


8 


-1.12721979 


-6.60059664 


107 


43 


5 


82 


237 


536 


9 


4.00046339 


1.92220990 


7 


36 


6 


51 


100 


301 


10 


4.22474504 


1.47400000 


37 


15 


195 


1173 


1420 


5344 


Average 














357.2 


943.2 



Table 11. Results of the experiments executed by the IBBA with the non-differentiable 
problems. 



Problem 


XIBBA 


FIBBA 


Ng, 




Ng, 


Nf 


Iterations 


Eval. 


1 


1.25830963 


4.17420017 


23 






28 


51 


79 


2 


1.95967593 


-0.07915191 


18 






16 


34 


50 


3 


9.40068508 


-4.40068508 


171 






19 


190 


209 


4 


0.33286804 


3.34619770 


136 


15 




84 


235 


418 


5 


0.86995142 


0.74167893 


168 


91 




24 


283 


422 


6 


3.76977185 


0.16666667 


16 


16 




597 


629 


1839 


7 


5.20120767 


0.90312158 


63 


18 




39 


120 


216 


8 


8.02835033 


4.05006890 


29 


11 


3 


21 


64 


144 


9 


0.95032461 


2.64804102 


8 


86 


57 


183 


334 


1083 


10 


0.79996352 


1.00023345 


42 


3 


17 


13 


75 


151 


Average 














201.5 


461.1 



- the column Nf shows how many times the objective function f{x) has 
been evaluated; 

- the column "Eval." is the total number of evaluations of the objective 
function and the constraints. This quantity is equal to: 

- Ng^ +2xNf,foT problems with one constraint; 
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Table III. Differentiable functions. Numerical results obtained by the PEN. 



Problem 


XPEN 


FXPEN 


P* 


Iterations 


Eval. 


i 


1 ACT 1 onn/1 
1 .Uj / 18UU4 


— /.DllOJOO/ 


1 c 
15 


OJ 


IDD 


L 


1 Ai ^r\mC/i 
i.UloUyzM 


D.4oi4Zoyo 


yo 


yj4 


lyuo 


3 


-5.99184997 


-2.94292577 


15 


119 


238 


4 


2.45953057 


2.84080890 


490 


1762 


5286 


5 


9.28468704 


-1.27484673 


15 


765 


2295 


6 


2.32334492 


-1.68307049 


15 


477 


1431 


7 


-0.77476915 


-0.33007412 


15 


917 


2751 


8 


-1.12719146 


-6.60059658 


15 


821 


3284 


9 


4.00042801 


1.92220821 


15 


262 


1048 


10 


4.22482084 


1.47400000 


15 


2019 


8076 


Average 








817.9 


2648.1 



Table IV. Non-Differentiable problems. Numerical results obtained by the 
PEN. 



Problem 


XPEN 


FXPEN 




Iterations 


Eval. 


1 


1.25810384 


4.17441502 


15 


247 


494 


2 


1.95953624 


-0.07902265 


15 


241 


482 


3 


9.40072023 


-4.40072023 


15 


797 


1594 


4 


0.33278550 


3.34620350 


15 


272 


819 


5 


0.86995489 


0.74168456 


20 


671 


2013 


6 


3.76944805 


0.16666667 


15 


909 


2727 


7 


5.20113260 


0.90351752 


15 


199 


597 


8 


8.02859874 


4.05157770 


15 


365 


1460 


9 


0.95019236 


2.64804101 


15 


1183 


4732 


10 


0.79988668 


1.00072517 


15 


135 


540 


Average 








501.9 


1545.8 



- A'gj + 2 X A'gj + 3 X Nf^ for problems with two constraints; 

- Ng^+2xNg2+3xNgj+4xNf, for problems with three constraints . 
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In Table ni (Differentiable problems) and Table IV (Non-Differentiable 
problems) the results obtained by the PEN are collected. The constrained 
problems were reduced to the unconstrained ones as follows 

fp* (x) = f{x) + P* max {gi (x) ,g2{x),. . . ,gN,{x),0} . (48) 
The coefficient P* has been computed by the rules: 

1. the coefficient P* has been chosen equal to 15 for all the problems and it 
has been checked if the found solution (XPEN,FXPEN) for each problem 
belongs or not to the feasible subregions; 

2. if it does not belong to the feasible subregions, the coefficient P* has 
been iteratively increased by 10 starting from 20 until a feasible solution 
has been found. Particularly, this means that a feasible solution has not 
been found in Table HI for the problem 2 when P* is equal to 80, for the 
problem 4 when P* is equal to 480, and in Table TV for the problem 5 
when P* is equal to 15. 

It must be noticed that in Tables 111, IV the meaning of the column "Eval." 
is different in comparison with Tables 1 and 11. In Tables 111, IV this column 
shows the total number of evaluations of the objective function f{x) and all 
the constraints. Thus, it is equal to 

{N„+\)X Niter, 

where A^v is the number of constraints and Nuer is the number of iterations for 
each problem. 

In Figures 4 and 5 we show the dynamic diagrams of the search executed 
by the IBBA and the PEN for the differentiable problem 7 from (Famularo, 
Sergeyev, and Pugliese, 2001): 

min fix) = exp (-cos (4x - 3)) + (4x -3)^-1 

xe[-3,2] 250 

subject to 

^ 1 

gi (x) = sin (x) exp(— sin(3;K:)) + - < 

g2{x) = cosQ(x + 3)^-sin(7(x + 3)) + ^<0 

The problem has two disjoint feasible subregions shown by two continuous 
bold lines and the global optimum x* is located at the point x* = —0.774575. 
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Table V. Differentiable problems: comparison between the IBBA and 
the PEN. 













Evaluations 


Prnhl r'TTi 


PEN 


IBBA 


V.) 1 LI LI L' 


PEN 


IBBA 


Speedup 


1 


83 


23 


3.61 


166 


36 


4.61 


2 


954 


227 


4.20 


1906 


248 


7.69 


3 


119 


62 


1.92 


238 


84 


2.83 


4 


1762 


953 


1.85 


5286 


1459 


3.62 


5 


765 


144 


5.31 


2295 


402 


5.71 


6 


477 


112 


4.26 


1431 


228 


6.28 


7 


917 


294 


3.12 


2751 


794 


3.46 


8 


821 


237 


3.46 


3284 


536 


6.13 


9 


262 


100 


2.62 


1048 


301 


3.48 


10 


2019 


1420 


1.42 


8076 


5344 


1.51 


Average 


817.9 


357.2 


2.29 


2648.1 


943.2 


2.81 



The first line (from up to down) of "+" located under the graph of the 
problem 7 in the upper subplot of Figure 4 represents the points where the first 
constraint has not been satisfied (number of iterations equal to 35). Thus, due 
to the decision rules of the IBBA, the second constraint has not been evaluated 
at these points. The second line of "+" represents the points where the first 
constraint has been satisfied but the second constraint has been not (number of 
iterations equal to 18). In these points both constraints have been evaluated 
but the objective function has been not. The last line represents the points 
where both the constraints have been satisfied (number of evaluations equal to 
241). The total number of evaluations is equal to 35 + 18 x 2 + 241 x 3 = 794. 
These evaluations have been executed during 35 + 18 + 241 = 294 iterations. 

The line of "+" located under the graph in the upper subplot of Figure 5 
represents the points where the function (48) has been evaluated. The number 
of iterations is equal to 917 and the number of evaluations is equal to 917 x 
3 = 2757. 

Finally, the infeasibility of the differentiable problem from (Famularo, 
Sergeyev, and Pugliese, 2001) has been determined by the IBBA in 38 itera- 
tions consisting of 9 evaluations of the first constraint and 29 evaluations of 
the first and second constraints (i.e., 67 evaluations in total). The infeasibility 
of the non-differentiable problem from (Famularo, Sergeyev, and Pugliese, 
2001) has been determined by the IBBA in 98 iterations consisting of 93 
evaluations of the first constraint and 5 evaluations of the first and second 
constraints (i.e., 103 evaluations in total). Naturally, the objective functions 
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Table VI. Non-differentiable problems: comparison between the the 
IBBA and the PEN. 







Iterations 






Evaluations 


Problem 


PEN 


IBBA 


Speedup 


PEN 


IBBA 


Speedup 


1 


247 


51 


4.84 


494 


79 


6.25 


2 


241 


34 


7.09 


482 


50 


9.64 


3 


797 


190 


4.19 


1594 


209 


7.63 


4 


272 


235 


1.16 


819 


418 


1.96 


5 


671 


283 


2.37 


2013 


422 


4.77 


6 


909 


629 


1.45 


2727 


1839 


1.48 


7 


199 


120 


1.66 


597 


216 


2.76 


8 


365 


64 


5.70 


1460 


144 


10.14 


9 


1183 


334 


3.54 


4732 


1083 


4.37 


10 


135 


75 


1.80 


540 


151 


3.58 


Average 


501.9 


201.5 


2.49 


1545.8 


461.1 


3.35 



were not evaluated at all in both cases. Note that experiments for infeasible 
problems have not been executed with the PEN because the penalty approach 
does not allow to the user to determine infeasibility of problems. 



6. Concluding remarks 

Lipschitz univariate constrained global optimization problems where both the 
objective function and constraints can be multiextremal have been considered 
in this paper. The constrained problem has been reduced to a discontinuous 
unconstrained problem by the index scheme. A Branch-and-Bound method 
for solving the reduced problem has been proposed. Convergence conditions 
of the new method have been established. 

The new algorithm works without usage of derivatives. It either deter- 
mines the infeasibility of the original problems or finds upper and lower 
bounds of the global solution. Note that it is able to work with problems 
where the objective function and/or constraints are not defined over the whole 
search region. It does not evaluate all the constraints during every iteration. 
The introduction of additional variables and/or parameters is not required. 

Extensive numerical results show quite a satisfactory performance of the 
new technique. The behaviour of the Index Branch-and-Bound method was 
compared to the method of Pijavskii combined with a penalty approach. This 
algorithm has been chosen for comparison because it used the same informa- 
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tion about the problem as the IBBA - the Lipschitz constants for the objective 
function and constraints. 

A priori the penalty approach combined with the method of Pijavskii 
seemed to be more attractive because it dealt only with one function. In the 
facts however, the evaluation of this function requires the evaluation of m + 1 
initial functions. The second disadvantage of the penalty approach is that it 
requires an accurate tuning of the penalty coefficient in contrast to the IBBA 
which works without any additional parameter. Finally, when the penalty ap- 
proach is used and a constraint g{x) is defined only over a subregion [c,d] 
of the search region [a,b], the problem of extending g{x) to the whole region 
[a,b] arises. In contrast, the IBBA does not have this difficulty because every 
constraint (and the objective function) is evaluated only within its region of 
definition. 
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Figure 1. Construction of the function (p(x) 
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Figure 2. The case v(x,_2) = v(x,_i ) = j, v{xi) = v(x,+i ) =7 + 1 
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Figures. The case v(;Ci_i) = y + 2, v{xi) = j 
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Figure 4. Optimization of the differentiable problem 7 by the IBBA. 
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Figure 5. Optimization of the differentiable problem 7 by the PEN. 
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